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ABSTRACT 

The kinetic theory of turbulent flow previously developed is 
enployed to study the mixing-limited combustion of hydrogen in axisymmetric 
jet. The integro-differential equations in two-spatial and three velocity 
coordinates describing the combustion are reduced to a set of hyperbolic 
partial differential equations in the two spatial coordinates by a 
bimodal approximation. The MacCormick's finite difference method is then 
enployed for solution. The flame length is longer than that predicted 
by the flame-sheet analysis, and is found to be in general agreement 
with a recent experimental result. Increase of either the turbulence 
energy or scale results in an enhancement of the conbustion rate and, 
hence, in a shorter flame length. Details of the numerical method as 
well as of the physical findings are discussed. 
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SUMARY 

Kinetic theory of turbulence developed previously was employed to 
analyze the turbulent diffusion flame in axisymnetric hydrogen jet. The 
fluid was assumed to be incompressible, and the mean flow field was 
considered to be given for convenience . The kinetic equations governing 
the chemical species and energy were transformed into a set of hyperbolic 
differential equations by the use of bimodal moment method which has been 
modified to circumvent certain' difficulty encountered in earlier works. 
These equations were integrated by a finite difference scheme. This 
method of solution was shewn to be flexible and applicable to many 
combustion problems of engineering interest. The solution was compared to 
the existing experimental data. Satisfactory agreements with the 
experimental data were obtained on the flame length and on the centerline 
variations of the chemical species concentrations. The maximum temperature 


of the flame was substantially below that predicted by the conventional 
flame-sheet solution. The overall combustion rate enhances and, 
therefore, the flame length shortens as the t ..<•>- julenee energy and scale 
were increased. Greater turbulence energy enhances the overall combustion 
rate by increasing both the mixing and the turbulent dissipation- 
limited chemical reaction rate. On the other hand, the larger scale 
gave a greater overall combustion rate through increased mixing, but a max- 
imum would be reached beyond which the combustion rate would decrease 
because of the reduced dissipation-limited reaction rate. 

INTRODUCTION 

Recognizing certain inherent limitations of the conventional 
method of analysis of chemically reacting flows, a kinetic theory of 
turbulence was developed previously. (See ref. 1-6.) Ibis theory 
describes the probability density function of the fluid elements 
comprising turbulence field. Probability density functions of the 
chemical species being transported by the fluid elements are also 
described in similar manner. 

The governing kinetic equations are lntegro-differential equations 
in both velocity and physical spaces, similar to those for the molecular 
kinetic theory, and present a formidable challenge to the available 
methods of solution. The only practical method of solution for chemically 
reacting flows was found to be the bimodal moment method, which was 
successfully employed to solve several chemically reacting flow problems 
of engineering interest (refs. 1-4). However, when the method was applied 
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to the combustion problems of initially unmixed reactants, it was found 
that the generated mean profiles contained certain discontinuities in 
the gradients which were not real. These discontinuities are direct 
consequences of the bimodal method of solution, as explained in ref. 5> 
and they do not affect the combustion rate and otter pertinent description 
of the physical problem. Nevertheless, a method should be found which 
will remedy this difficulty. 

The ultimate answer lies in the exact numerical integration of the 
kinetic equations. The equations were solved exactly at least in one 
velocity space for a combustion problem in ref. 7, and the solutions gave 
correct continuous mean profiles to all orders as expected. The numerical 
complexity, however, was such that the solution was confined to a 
relatively simple shear flow. 

As an interim remedy of the difficulty with blmodal method, an 
additional kinetic equation was constructed in ref. 6 which governs the 
correlation function representing the chemical reaction. This equation 
circumvented the difficulty of the blmodal moment method, and the method 
was successfully applied to solution of the turbulent chemical laser 
problem in ref. 6. 

Aside from the problems discussed above, there exists the difficulty 
of numerically Integrating a set of hyperbolic partial differential 
equations resulting from the moment method in two physical dimensions. The 
MacCorroick’s finite difference scheme (ref. 8) was shown to be useful in 
the solution for simple two-dimensional flow (ref. 6). However, application 
of this scheme to the axisynmetric jet turned out to be a formidable task. 
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A substantial portion of the effort expended In the present work 
concerns the numerical Integration of the governing moment equations 
by the MacComick ' s method. 

Text of this report begins with formulation of the kinetic 
equations governing the mixing-limited combustion in axlsyimetric jet. 
The fluid was assumed to be incompressible and the mean flow field 
was considered to be known in terms of the existing incompressible 
solutions. Appropriate derivation and transformation of the moment 
equations are given. Solutions for the hydrogen combustion in air 
are then obtained and compared to the available experimental data. 
Because of the importance of the numerical scheme to the success of the 
kinetic-theory approach to turbulence, a detailed description of the 
numerical solution technique is included in the appendix. 

SYMBOLS 


a stoichiometric coefficient for oxidant 

a^a^a^a^ constants for Eqs. (14), (16), and (17). 


b 

c 

c i 

C P 

D 

d 

E 


stoichiometric coefficient for fuel 
mass fraction 

i-th mass fraction used in appendix 
specific heat 


molecular diffusivlty 


stoichiometric coefficient for combustion product 
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P function defined in Eq. CA^J 

f probability density function of fluid elements 
h C p T/Ah° 

Ah 0 heat of combustion per unit mass of product 
kj, forward reaction rate coefficient 

L nozzle diameter 

M molecular weight 

Q general function of 9 

R r/L 

r radial coordinate defined in Fig. 1 

S function defined in Eq. (A5) 

T absolute temperature 

U x-c exponent 0 f instantaneous velocity relative to the mean velocity 

i-ccmponent of instantaneous velocity relative to the mean velocity 

U instantaneous velocity vector relative to the mean velocity vector 

u x-ccmponent of instantaneous velocity 

u. average center line velocity 

u i i-component of instantaneous velocity 

Up average velocity at the nozzle 

u w average free stream velocity 

V r-component of Instantaneous relative velocity 

v r-component of instantaneous velocity 

W instantaneous relative velocity component normal to x-r plane 
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Instantaneous velocity component normal to x-r plane 
x/L 

axial coordinate defined in Pig. 1 

general mass fraction representing c ft , c b , etc. 

functions defined by Eqs. (22) - (24) 

functions defined by Eq. (37) 

« / fQ dU 

dissipation ratio in turbulence 

1 - <vv 

eddy diffusivity 

similarity variable defined by Eq. (17) 
turbulence scale 

X/L 

characteristic directions defined by Eq. (39) 
functions defined by Eqs.(A9) and (A10) 
density 

constant for Eq. (17) 
function defined by Eq.(32a) 

functions defined by Eq. (36) 
functions defined by Eq. (32b) 
functions defined by Eqs. (A3) and (A6) 

instantaneous production rate per unit mss by chemical reaction 



Subscripts 


a oxidant 

b fuel 

c center line 

d product 

e chemically inert species 

h °h s h 

i,j ,k Cartesian tensor indices 

o ensemble averaged 

p at the nozzle 

00 free stream 

Superscript 

vectors 

Note that all symbols for velocities, except Up , u, , and u w , 
represent dimensionless quantities after Eq. (26). 

FORMULATION 


Chemical Reaction 

The standard axisyirmetric turbulent jet shown in Fig. 1 is studied. 

The Jet, consisting of a fuel and a chemically inert species, is surrounded 
by air, and the following one-step combustion is considered to take place. 

k f 

a(oxidant) + b(fuel) ► d(product) (1) 
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Prom the law of mass action, the rates of production by chemical 
reaction of the three reactants are given by 


M - ir d ” d « a+t> -l „a b 
w d “ k ftfT3> p a b 



( 2 ) 


All symbols are defined in the nomenclature. 


Kinetic Equations 

Starting point of formulation of the governing equations 1 b the 
kinetic equation (see refs. 1, 5, and 6), 


(u, 


Oj 


+ Uj ) 


3fZ 


<u kV 


<w 

— 2X 


1/2 


(1 + y) 


5UJ (fU J z) 


3 2 fz 

* 5*5 


]- 


<u k u k > 1/2 Y 

2A 


f(z - z 0 ) 


+ fu>_ 


(3) 


where j and k are the C&rtesian tensor indices. The symbol z stands 
for difinlte quantity being transported by the fluid elements. For the 
combustion described byEqs. (2), Eq. (3) represents the five kinetic 
equations respectively for 
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z 


°a» V c d» c e 81x1 h 


W 


These equations will be solved by the blmodal moment method 
subsequently. As was discussed In the proceeding section, certain 
discontinuities in the gradients of the mean profiles are expected 
in the results when ■* ». Although these do not affect the combustion 
rate and other physical descriptions of the problem, the discontinuities 
themselves are unreal consequences of the blmodal method of solution. 

As was discussed In ref. 5, the blmodal method approximates the 
continuous distribution of a reactant concentration among the fluid 
elements with various instantaneous velocities by two discrete values of 
the concentration depending only on t,.e sign of V. At the edge of a 
diffusion flame zone, therefore, concentration of one of the reactants 
contained in all fluid elements becomes zero simultaneously. This results 
in the discontinuous gradients of the mean reactant concentrations at 
the flame edge. Such unreal discontinuities naturally disappear as more 
exact method is employed for solution of the kinetic equations (refs. 5* 
and 7). These methods (refs. 5 and 7), however , are of such complexity 
that their application to practical flow problems such as the axisynmetric 
Jet is extremely difficult. 

As an interim technique to circumvent the difficulty with the blmodal 
method, until a acre viable method of solution is found, it was suggested 
in ref. 6 that an additional kinetic equation be employed. This 
equation is constructed by letting 




a b 
2 * e a°b 


(5) 


In Eq, (3) as 


3fcjc5 


<W 


1/2 




«W ^*44 1 -,.l b a .b . 

”1 — 1515- J ?X f(c a c b “ c ao c bo ) 


* f«ab* 


( 6 ) 


Expression for is constructed such that the first moment of 

Eq. (6) corresponds to the conservation equation for < c®c^ > derived 
fran the Navler-Stokes and standard species conservation equations. When 
a ■ b ■ 1, the expression for ./as obtained in ref. 6 as 

w ab " °b w a + c a w b (7) 

such that the first moment of Eq. (6) was term-wise comparable to the 
conservation equation for <c & e b > derived fran the standard species 
conservation equations for c & and c^, subject to the correspondence 


3 c BCj, 

A — s 

»v 




2X (<c aV " c ao c tio )l 


( 8 ) 


The presort analysis is mainly concerned with the combustion of 
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ijOki.... 


hydrogen 


0 2 ♦ 2Hg 


2*^0 


(9) 


for which a ■ 1, and b • 2. For these stoichiometric coefficients, 
it was found that • * 


w ab " °b w a + 2c a c b w b (10) 

gives the tem-wise compatibility of the first moment of Bq. (6) with 

2 

the conservation equation for <c & c£> derived from the conventional 
conservation equations for c ft and c fc , subject to the following 
correspondence for dissipation, 


^ ^ 5^ > * < C «55J 3JJ>) 

■ ^ M-vS) 


(id 


The dissipation correspondence of Eqs. (8) and (11) are plausible 
hypotheses, and consistency of Eq. (6) with the original kinetic equations 
for e , c^, etc. cannot be established. Therefore, the bimodal 
solutions using this additional equation must be considered interim 
ones until more exact solutions of the integro-differential kinetic 
equations can be obtained. 
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MixlngyLlmlted Reaction 

When the chemical reaction rate, k f , is sufficiently larger 
than the mixing and other characteristic rates of the problem, Eqs. 
(2) demand that 



0 


( 12 ) 


for all x and U. Eacl term of Eq. (6), then, vanishes except the 
last two terms, the dissipation and chemical reaction terms, which 

become with the use of Eq. (10) and Eqs. (2), 


foi, 


a 


<U k U k > 1/2 Y 

— 2S 


c c 2 
c ao c bo 


(4 


■V 

FT c 


a 


a c b. 


(13) 


As explained in ref, 1 , the above shows that, in the limit of fast 
reaction, the reaction rate is equal to the turbulent dissipation rate. 

The set of kinetic equations governing the hydrogen combustion 
now consists of the five equations represented by Eqs. (3) and (4), and 
the degenerated equation, Eq. (13). Solution of this set is discussed 
in the next section following the brief consideration of the flow 
field given below. 


Flow Field 

Mean flow field for incompressible turbulent Jet is well known. 

(See refs. 9 and 10.) So that full attention can be given to describing 
the combustion phenonraenen, the mean flow field solutions found in 
refs. 9 and 10 are employed in the present analysis instead of solving 
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for the flow field by the kinetic theory. 

Recognizing the fact that exact details of the flow field have 
only a secondary effect on combustion, solution Of ref. 9 for u^u^ = 0 
is first modified slightly to include the influence of the finite initial 
jet diameter on the flow field using the information given i ref. 10. 

The result is then generalized so that it can be approximately used for 
arbitrary values of u also. The mean velocity components employed 


where 
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c 
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e - 1 
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= 6.4, a 2 55 0.6, and a^ - 
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Having completed forma lat ion of the problem, solution of Eq. (3), 


with Eq. (4) , and Eq. (13) constitutes the next section. 


ANALYSIS 


As stated earlier, the kinetic equations governing combustion will 
be solved by the bimodal moment method. The method is similar to that 
employed in the previous work, ref 6, except the mean velocity component, 
v Q , is not zero in the present problem, and this requires a somewhat 
different approach in the division of the velocity space. 


Moment Equations 


Ihe first two moments of Eq. (3) are obtained by multiplying the 
equation by 1 and V successively and by integrating the resulting 
equations term-wise with respect to the velocity space. In cylindrical 
coordinates, these equations are 


+ u)f 2 ® 


V + rh r f (v o + v)f z = Jf w z dU , 


u o h / V 


Vf z d U + v 


o 3r 


Vf z dU + 


I ?/ Vt 2 


/ 3v 3v A r + 

( u o 3 r + v o wjj f 2 ® 


+ r f 2dU 


<u k u k >V2 C - f 

— (1 + 2y) J Vf z dU +J Vf u> z dU 


ifl 


In the derivation of Eq. (19), the continuity equation resulting from 
setting z - 1 in Eq. (18), and the relationships 


I UVf z dU < < uj Vf z dU , 


Uf z dU < < 


u 0 Jf z dU , 


have been employed. 


Bimodal Approximation 


For convenience , the probability density function of the fluid elements 
is approximated by a Maxwellian function in terms of the velocity components 


relative to the mean velocity. Thus, 


U 2 + V 2 + W 2 


* <U kV> 


The species mass fraction, z , is then approximated by a bimodal 
function as. 


z(x,U) = z(x,r,V) 


= z,(x,r) + z 0 (x,r), 


Mltta 


MMMp 


Note that the above relationships Imply 


z 2 (x,r) * 0 for v > 0 , 

and ( 

z^Xjr) «* 0 for v < 0 . 


Therefore, the bimodal division is actually carried out in the absolute 
velocity space. 

Eqs. (21) and (22) facilitate the evaluation of a moment as 


oo oo 


<Qz> 


= y*Qf z dU = j dV J dW J dU Qf 


— V —oo 

O 




dU Qf z , 


—00 —00 


where Q denotes an arbitrary function of U. . 


Governing Equations for Combustion 

With the use of Eqs. (21), (25), and (13), and the relationship 
between , and u) & given in Eqs. (2), Eqs. (18) and (19) are 

reduced to a set of explicit partial differential equations in x and 
r for each of the z defined by Eq. (4) . Then, with appropriate non- 
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dlmenaionalization and some manipulation, the following four sets of 
equations are ootained. 

In all subsequent equations, the symbols for velocities, except 
, u c , and , represent the normalized velocities with respect 


u 0 " V u c * 

T o ' V u c ■ 




U = U/u c , 

V • V/u 0 , 

and 

W = W/u . 
c 

(26) 


« c a > 


f(c a ) - 


E^v c ao°bo 
2\ 2 


2 {♦+ + f erf v 0 j | 

) c ao + <V V S^] 


i 3v u 3v rt 

_2. 0 c + v 2. i c + <vc > 

a 3X + v o 3R ) c ao + vc a 
c 


(27a) 


? V2 


2 r (1 + 2Y)<v ^ 


exp 


(- A v 2 ) 

\ 2E v o / * 


( 27 b) 


*(%) = - jf 


^ E 1/2 y c ao c bo 


M a X 2 


[* + + *' err [(si ) 12 v o]| - 


(28a) 


nc b ) 


r/u 

3v u 

3Vq \ 

3v 1 


o c 
3X 

+ v o 3R"j c bo + 

<Vc b > 3T J 
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gl/2 Ml r? 

IT (1 + 2lr)<v V - *1 


c ao°bo *' exp 


(■ 5 § v o ) 


M, pl/2 c c 2 ( , 

a/ ^ n . u E Y SO DO J a t ± x m 

*(c d ) - H s-i* 5 — < $ + ♦ 

a a \ 


*(«d> 


K u. 3v_u. 3v \ 

tg- hr '\ ifr ) 


c ao + 
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c ao c bo exp 
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+ 4> — erf 
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The functions governing the chemically Inert species are obtained from 
those of the other chemical species as 


f c el dU « f 


: al + c bl + c di : 


f c e2 dU 


IU » f f 1 " (°a 2 + c b2 + c d2^] dU * 


The operators $ and 4* are defined by 


1(z) ■ u o ir +v o?t + rIk' r<Vz> > • 


,(,) ’ t 


u 3u <Vz> 
o c 


3<Vz> . 3 

'oir + 3R 


L <V 2 z> . 


Various ensemble averaged quantities appearing in Eqs. (27) - (32) are given 


1 1 

z Q * <z> » g ( z i + ^ 2 ) + j (z 1 - z 2 ) erf 


KC ■•] • 


! 1 ' Z 2 } (sf) 


6XP 


<V 2 z> » | (z, +z~) + (z, - z-) j | erf 



(35) 


-(^) 1 / 2 v o eItp (-^ v o) j 


Functions 4> and <f arising from Eq. (13) are defined as 


fa) + °b y + , + -w + 

\ 2 — / 5^ (c a + c a )(c b + c b ) 




r <0 a - c a Ko b - °b> 


(36) 


where 


■ z x + z 2 , 


and z" 


Z l- Z 2 


(37) 


With the velocity profiles given by Eqs. (14) - (17), and the 

A 

turbulence properties, E and X , considered to be known in the present 
analysis, as will be discussed subsequently, Eqs. (27a) through (3D constitute 
ten equations for the equal number of unknown functions, c^, c^, c^, c^, 
c dl* c d2* c e^* c e2* ^1* ant * ^2 ’* Solution of these equations and discussion 
of the results comprise the remainder of this report. 


•One may alternatively consider c*, c a , c£, c^, etc., as the 
unknown functions via Eqs. (37) - 
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Wave Behavior of Equations and 


SEE 


Conditions 


As was stated earlier, the governing equations, Eqs. (27a) - (31), are 
hyperbolic. These equations will be solved by a finite difference scheme j 
however, it will be useful before the solution to study their behavior on 
the characteristic plane. Such a study, among other things, will clarify 
the boundary conditions to be applied. 

The sets of two equations, (27a) and (27b), (28a) and (28b), etc., are 
all similar. Equations (27a) and (27b), therefore, are used in the following 
discussion of the general behavior of the governing equations. 

The terms <{f erf v Q in Eq. (27a) and [ (u o /u c )(3v o u c /3X) + 

W 3R) ] c ao 111 Eq ‘ (2Tb) do not affect the basic behavior of the equations . 
Neglecting these terms for the present discussion, Eqs. (27a) and (27b) can 
be transformed into the following set of characteristic equations. 




/2EN 172 1 - T f 1 

W R a L8M- 


+ UJ c a u o dX 


E^ 2 y +, + 






H u o wJ 


For the axisynmetric jet shown on Fig. 1, the characteristic 

directions defined by Eq. (39) are sketched In Fig. 2. The two waves 

+ 1/2 — + 1 /2 — 
consisting of [c_ - (2/ir) cj and [c + (2/ir) c„] propagate along 

the two characteristic directions £ + and £” respectively according to 

Bq. (38). Ihese waves dissipate by the quantities contained in the last 

curly bracket of Eq. (38). These quantities represent three physical 

phenomena. The first and third terms involving <J> + and f denote 

the chemical reaction, the second and the last terms signify the turbulent 

mixing, and the term next to last represent turbulent dissipation. These 

physical phenomena combine to dissipate the waves as they travel along the 

characteristics. The wave behavior will be discussed further following the 

numerical solution later. 

As seen from Eq. (39) and Fig. 2, the R-axis is "space-like" whereas 
the X-axis is "time-like" (see ref. 11). Therefore, two boundary conditions 
are required along the F-axls and one along the x-axis. These requirements 
can be satisfied by the conditions that, 

for x ■ 0 and all R * 0 , 


a 


^x-0 , 


(*K>) 


and 


a 


^a^x-0 * 


m 


For R - 0 and all x ^ 0 , 
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c“ - 0 . (42) 

By the use of Eq. (33) » Eq. (40) can be replaced by the mean mass fraction 
profile, boundary condition of Eq. (42) arises from the 

reflective wave condition, 



which can be satisfied only by Eq. (42). 

SOLUTION 

As was stated previously, the governing equations, Eqs. (27a) through 
(31), are solved by the MacCormick’s finite difference method (ref. 8). Details 

of the numerical method and the computer program are given in the Appendix . 

The turbulence properties and the boundary conditions corresponding to Eqs. 

(40) - (42) employed are discussed in this section. 

Turbulence Properties 

Turbulence energy and scale are needed to render the governing equations 
self-contained . Here, these quantities are chosen such that they will 
give the conventional eddy diffusivity (ref. 9) in the limit of the 
statistically near-equilibrium turbulence field (refs, 5» 6). As was done 
for discussion of the wave properties of the governing equations in the 
preceeding section, Eqs. (27a) and (27b), with the term l(u /u )(3v u_/3X) 

O C O v 
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+ v 0 Ov < /*R)lc ao neglected, are used for demonstivticn, General 
behavior of all sets of the two equations, (27a) and (27b), (28a) and 


(28b), etc., is the sane. 

Equilibrium turbulence field is approached as 


X 


m 


In this limit, Eqs. (27b), (3*0, and (36) show that 


<Vc a > 

a 


a 


and 


0 , 

0 , 

0 . 


( 45 ) 


Eq. (27b) then reduces to the conventional , gradient driven diffusion 
relationship. 


<Vc > ■ 
a 


2 7 r.1/2 3c ao 

Tnrrm x E ir 


(46) 


The eddy dlffusivlty is defined by the above equation as 



KtHyT 




(47) 


A substitution of Eq. (46) Into Eq. (<^7a) results in the conventional, 
parabolic conservation equation (refs. 9, 10) with the same flame- sheet reaction 
(ref. 12), 





1 

r 


3 

IF 


( r S,Tr) + “«o S(r - r *> 


(«) 


where r* denotes the position of flame sheet. 

The turbulence scale and energy, \ and E , are chosen to satisfy 
Eq. (ii/) with the conventional value of eddy diffUsivity, e m . The mean 
velocity profiles to be erployed, Eqs. (14) - (17), have been obtained from 
ref. 9 for the axisyirmetric jet with u„/Up ■ 0 . The same reference gives 

o - 15.17 , 

and (49) 

e ro - 0.01388 L . 

Implicit in the conventional similarity analysis of ref. 9 is that all length 
scales vary as l/(a ^ * X) according to the variable rj given in Eqs. (17). 

A 

Therefor# t the present turbulence scale A should be a constant fraction of 
the Jet radius, and should vary with respect to X such that n is a constant 
Denoting this constant by , Eqs. (17) give the relationship 


\ 



( 50 ) 


A substitution of Eqs. (14), (49), and (50) into Eq. (47) results in 


0.0493 (1 + 2y) 






9 


The turbulence dissipation ratio, y , is uet to 1/3 as was done In 
the previous work, ref. 6 . 

Eg. (51) shows that E must be a constant valr-' in order to produce 
the conventional conservation set, Eqs. (47) and (48), in the limit of near- 
equilibrium turbulence field. Existing literature (see, for instance, refs. 10 
and 13) shews that E is of order 0.1. The values of E employed in 
the present computation are, 

E ■ 0.05 to 0.2. (52) 

The corresponding values of a^ are determined from Eg. (51). Note that 
ajj ■ 0.26 when E ■ 0.1. Since the Jet boundary corresponds to n % 2.5 , 
the turbulence scale represented by a^ ■ 0.26 is about one tenth of the jet 
radius which is of the order of the mixing length as expected. Solutions for 
E ■ 0.1 and a^ ■ 0.26 will be used for comparison with the experimental 
data subsequently. 

Boundary Conditions 

The boundary conditions for Eqs. (27a) - (30b) are the same as Eqs. (40) 

- (42) for each of the variables, c* , c“ , c£ , , etc. The conditions 

of Eqs. (40) and (41), however, also imply that each of the variables is 
specified at R • for all X . 

Before the boundary conditions are explicitly specified, a consideration 
should be given to the following. The similarity velocity prattles given by 
Eqs. (14) - (16) are valid for X £ 5*8 (see ref. 10). In fact, ususally, 
the flow is not fully turbulent for X < 5.8 and there the present turbulence 
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m 
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analysis does not apply. Since the present interest lies with the 
combustion description for large X’s where influence of the initial 
condition is negligible, the initial conditions are specified at X » 5.8 . 



Solutions of Eqs. (27a) - (30b) are obtained for the following 
initial and boundary conditions . 


5.8 , 



c ao “ c ao (R) * 

\ . 

(53) 

°bo “ c bo (R) » 

% " • 

(54) 

A 

o 

II 

•8 

o 

oP 1 

II 

o 

(55) 

h 0 - h 0 (R) , 

h“ * h“(R) . 

(56) 


For R -*■ ® and X > 5.8 , 





The Initial profiles for and c bQ shown on Fig. 3 are specified 

for Eqs. (53) and (54). The profile for c bQ is chosen such that it 
corresponds to the pure hydrogen Jet issuing from a circular nozzle of 
radius 1/2 as seen in Fig. 3; that is 




( 62 ) 



The c ac profile of Fig. 3 is chosen arbitrarily with oily two conditions 
in mind. The mixing with hydrogen begins at X = 5.8 , and that c profile 

dO 

is continuous with Eqs. (57) satisfied. Any overt discontinuities will tend 
to propagate along the characteristics and may cuase numerical difficulty. 

An arbitrarily continuous profield for h Q is specified at X = 5.8 for 
Eqs. (56) to satisfy Eqs. (60). 

The simplest value for c~ , c b , and h” to be assigned at X = 5*8 
is zero. It was found, however, that numerical stability of the solution is 
improved by specifying a small but non-vanishing function for each to satisfy 
Eqs. (57), (58), and (60). Both c dQ and c^ are specified to be zero for j 

all R initially. 

As was stated earlier, the present interest lies with the combustion 
description for large values of X and with the overall combustion phenomena 
for the entire jet. Details of the initial conditions have very little 
effect on these. 

Solutions of the governing equations, (27a) through (31), satisfying the 
boundary conditions, Eqs. (53) - (6l), are obtained by the McCormick's (ref. 

8) finitie difference method. Details of application of the method are found 
in the Appendix . Results of the solution are discussed in the next section. 
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DISCUSSION OP RESULTS 


Pigs. 4, 5j and 6 shew the mean reactant mass fraction and temperature 
profiles for the three different axial locations of X = 21.81 , 61.82 , and 
105.81 respectively. These are the solutions for E = .1 and a^ = 0.26 
which together comprise the eddy diffusivity given in ref. 9 in the limit 
of an equilibrium turbulence field (see the preceding section) . Therefore, these 
solutions are considered as the more representative ones for the standard 
jet issuing into a stationary air. 

Solutions were found to be insensitive to u^u^ when it is smaller 
than 0.2. On the other hand, the numerical difficulty increased rapidly 
with decreasing u^/Up . Therefore, the solutions discussed in this section 
were obtained for = 0.2 . 

The general shapes of the mean profiles are same as those found in 
the experimental works of refs. 14 - 17 . Note that the present results, 
c ao and Cp 0 , are in mass fractions whereas those given in refs. 14 and 15 
are in mole fractions. Mean profiles of the combustion product are similar 
to those of temperature , and are not shown in the figures. 

The temperature peak is located at a point away from the jet axis for 
small values of X (see Fig. 4). The peak point is slightly on the fuel-rich 
side from the stoichiometric location, where the fuel-to-oxygen mass ratio 
is 1/8 , in agreement with the finding of ref. 14. Location of the temperature 
peak moves to the jet axis as X is increased. 

At the larger X's (Figs. 5 and 6 ), oxygen survives through the 
combustion zone and appears in the jet axis. The flame comprising a shell 
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in the Jet (Pig. 4) for smaller X’s pervades throughout the Inner core 
of the jet (Pigs. 5 and 6) for larger values of X . The flame thickness 
is of the order of one-tenth of the jet axial distance, X , and comprises 
about one-half of the jet cross section. 

Distribution with respect to R of the combustion rate is shown in 
Pig. 7. Note that the maximum burning point does not correspond to the peak 
temperature position; particularly for X = 61.82 and 105.81 . For all 
X's , the maximum combustion point is slightly on the fuel-rich side of the 
stoichiometric position. For X = 21.81 , the maximum burning takes place 
at a point between the temperature peak and the point of stoichiometry. 

Comparison of Pigs. 8, 5, and 9 shows the effects of varying 
turbulence energy on the mean profiles. Other than the substantial change 
in the burning rate evident from the level of c bQ , E does not basically 
change the mean profiles. The combustion rates will be elaborated upon 
subsequently in conjunction with another figure. Thickness of the flame 
is seen to be insensitive to the turbulence energy variation. 

Figs. 10, 5, and 11 show that the flame thickness is substantially 
affected by the turbulence scale represented by a^ (see Eq. (50). Greater 
the scale (size of the energy containing eddies), greater is the flame 
thickness. This was expected from the previous studies (refs. 1 and 6). 

The flame is thick in contrast to the infinitesimally thin flame sheet 
predicted by the conventional theory (see ref. 12). Therefore, the flame 
shape and length for the jet flow cannot be precisely defined. However, 
loci of the stochiometric locations may be used for comparison of the 
flame shapes and lengths. Fig. 12 shows the stoichiometric positions 
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for the four of the computed cases. Both Increased turbulence energy and 
scale result in the shorter flame lengths. 

Shorter flame lengths, of course, are caused by the enhanced burning 
rates, and this is seen in Pigs. 13 and 14. These figures show the total 
integrated burning rates across the jet cross-section as functions of X . 
Stoichiometric lengths of the flames for E ■ 0.1 and 0.2 with a^ = 0.26 
found in Fig. 12 are also shown in Pig. 13. Because of the very thick nature 
of the flame zone, the combustion rate does not become zero at the 
stoichiometric end of the flow (Pig. 13), but it asymptotically approaches 


zero. 


Both larger turbulence energy and scale result in enhanced combustion 
rates (Fig^. 13 and 14), and, hence, in the shorter flame lengths (Fig. 12). 
Increased turbulence energy enhances both mixing and the chemical reaction 
rate which is controlled by the turbulent dissipation (Eq. 13). This is 
evident in the characteristic equations, Eqs. (38) and (39). Larger 
turbulence energy implies, according to Eq. (39), the greater slopes of 
the characteristics which cause a faster mixing between the fuel and the 
surrounding air. The increased turbulence energy also enlarges the 
dissipation-limited chemical reaction terms of Eq. ( 38 ), which are the 
terms involving <j> + and . Therefore, turbulence energy has a 
substantial influence on the combustion rate and flame length. This is in 
agreement with the previous findings (ref. 6) in plane mixing layer of infinite 
domain. It was found in ref. 6 that the combustion rate increased with . 


Although the present jet problem is quite different from the plane mixing layer 
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problem analyzed in the reference, it is interesting to note that the two 
Jet combustion rates respectively for E ■ 0.1 and 0.2 , as represented by 
reciprocals of the two stoichiometric lengths of the flames (see Pig. 12), 
scale very closely as vl? . 

Larger turbulence scale, X , enhances the combustion rate through 
increased mixing which is brought about in a manner quite different frctn 
that for the larger E discussed above. The characteristic directions 
are Independent of the turbulence scale (Eq. (39)). Increased scale 
actually reduces the turbulence dissipation- limited chemical reaction 
rate as seen in Eq. (38). In fact, magnitudes of all the terms on the 

A 

right-hand side of Eq. (38) are reduced by the Increased X except the 
mixing terns which are the second and last terms in the second curly 
bracket. However, as explained in the sentences following Eqs. (38) and 
(39), this ri^it-hand side constitutes the wave dissipation.* With the 
dissipation action on the wave diminished, the wave propagates further and 
for a longer period of time. For instance, the wave consisting of the fuel 


properties propagates into the oxygen-rich territory until the chemical 

reaction, the dissipative action on the wave, completes its task. Therefore, 

pi a 

. increased X simply stretches the reaction zone as seen in Figs. 10, 5, and 

I i 11, and in so doing, brings the reaction zone into the oxygen-rich region. 



•Note that the word "wave dissipation" is employed to mean the standard 
dissipative action on the wave. This should not be confused with the 
"turbulence dissipation" which refers to the viscous dissipation of 
the turbulent fluctuations. 
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This implies an enhanced mixing and, consequently, an enhanced combustion 

i. 

rate. This is in agreement with the findings of ref. 6 for the plane 
mixing layer.. As was found in that reference, however , it is expected that 
continuous increase of the scale will eventually result in the reduced 
combustion rate. This is because a point must be reached with increasing 

A 

A beyond which no amount of efficient mixing can overcome the inefficiency 
of the turbulence dissipation-limited reaction rate caused by the large 
scale. 

Typical variations of the reactant mass fractions and temperature 
along the Jet centerline are shown in Figs. 15-17. As explained earlier, 
the temperature peak occurs away from the centerline for small values of X . 
The maximum temperature is also indicated in the figures when it does not 
correspond to the centerline value. When the conventional flame-sheet 
solution is obtained, the temperature peak exists at the sheet which is located 
away from the centerline until the end of the flame. Also, this temperature 
is constant until the flame ends and it can be readily computed . This flame 
sheet terrperature is shown in the figures for comparison. Because the 
chemical reaction is spread across a thick flame zone, the peak temperature 
computed in the present study is substantially below the hypothetical 
flame sheet temperature . This difference is greater than the difference 
found in the experimental work of ref. 14. As was stated previously, 
variation the combustion product is similar to that of the temperature. 
Since concentration of the combustion product was measured in ref. 15 , but 
not the temperature, as well as in ref. 14 for pure hydrogen jet combustion, 
the present computed values of the combustion product is compared with those 
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of refs, 14 and 15 in Fig. 18 . The present values of the HgO concentration 
agree rather closely with the experimental values of Hawthorne, et al, (ref. 15) 
lift are considerably below those of Kent (ref. It). 

Fig. 18 shows that the computed concentrations of Kg and 0 2 agree 
satisfactorily with the experimental data of both refs. 14 and 15. Note that 
the 0 2 centerline concentration was found to be practically zero in ref. 15 
up to X of about 66 , which is in agreement with the present results. 

Finally, the computed stoichiometric contour of the flame is compared 
with that measured in ref. 14 in Fig. 19. No detailed agreement in these 
contours was expected since the present solution was based on an 
incompressible flow assumption. However, it is rather surprising to see 
the close agreement on the flame length, that is, on the overall combustion 
rate. Hawthorne , et al, (ref. 15) measured the "luminous" flame lengths 
for various jet flow conditions and found them to be between about X £ 128 
-*■ 146. These values also agree closely with the computed results. The 
computed contour for the hydrogen mass fraction of 0.003 is also shown 
in Fig. 19. This shows that the flame is thick and the combustion continues 
for a substantial distance following the stoichiometric end of the flame, though 
at diminishing rates, as explained earlier. 


CONCLUDING REMARKS 


The kinetic theory of turbulent flow was employed to analyze the mixing- 
limited combustion in axisymmetric jet. The bimodal moment method of 
solution of the kinetic equations was modified to circumvent certain 
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difficulty encountered previously. In the absence of a more direct 
numerical capability for solution of the original kinetic equations, this 
modification enables one to obtain satisfactory solutions of the turbulent 
combustion problems. It is noted that a direct numerical technique is 
being tried elsewhere (see ref. 7) with some success, and such a method 
should eventually replace the present method. 

The modified blmodal method was employed to transform the kinetic 
equations governing the combustion energetics into a set of hyperbolic 
differential equations. Assuming, for convenience, that the incompressible 
mean flew field is given, this set of the equations was integrated by the 
MacCormlck ’ s finite difference scheme for combustion of hydrogen. 

Detailed structure of the flame as well as the overall confcustion 
rate of the hydrogen jet were computed. The centerline variations of 
the reactant concentrations and the flame length obtained were compared 
with the available experimental data which showed a satisfactory agreement. 

A full discussion is found In the proceeding section. 

Che of the contributions of this report is the successful application 
of the finite difference scheme to solution of the governing equations for 
combustion based on the kinetic theory. This numerical technique — together 
with the modified moment method — should provide a considerable flexibility 
for analysis of the various combustion problems of engineering interest. 
Details of the application of the numerical method are given in the appendix. 
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APPETfXtX 


INTEGRATION OP THE GOVERNING EQUATIONS 


Eqs. (27a) - (30b) are transformed into a first order vector 
differential equation. Numerical solutions are obtained for the variables 
cj and c" where 

+ _ + + + + + + + 

C 1 " c a » c 3 " c b * c 5 " c d » c 7 " °h • 


c 2 " c a » c 4 " °b * c 6 " c d * c 8 " °h 


The vector differential equation is 


|| + -g|p (B) - S (fi) . 


Conponents of the vectors are given as follows. 


For i ■ odd; 


% ■ °i • 


F 1 (S) ■ ^ °io + 4 <Vo i > - * ^ exp (s$ v o ) °io • 
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LV^ar-^srVO 


+ Br] «p ( J V* )<vcj> 

+ Jr** (£*!)£ . tt 

where 

<! - m* hM .1* Of «,. (27.), (28a). (29.), or (30.) (« 

ana w 

* ■ Right M .I* of B,. (J7b)> (28b)> (J9b)> ^ 

3v„ 

• BT <V V • 

(A10 

undlrferontlated f. S(fi) 1, treated „ Bm( ^ 

numerical method employed 4 «a v,, Ba j 

™vj.oyea is based on the two-st*n a^n.u 
k r two-step explicit procedure 

developed by Lax and Wendroff (ref ift ) 9rv o 

‘ 18 and * more Gently, by ffecCormick f~r 
8)- 'The intermediate state n. < (ref. 

state, , is computed as 


cn the two-step explicit procedure 


computed as 


°J - ffi [ f <^ +1 ) - W^) ] + ax S(B) . 


flnal state 18 then determined by 

"5* 1 ’ 0 -*{? + v-w[f(S 


I,*) - W, ) 



♦ AX S(flj«)| . 


(A12) 


- -4 


Thus, for a given vector ? at X* nAX , the subsequent state vector is 
,'nr-nuted forward in the X-direction. The two-step explicit procedure provides 
a solution with second order accuracy in the finite difference increment. This 

i 

method has an advantage in which no mextrix inversion and iteration are 
required. It has, however, the disadvantage of being conditionally stable. 

Stability condition of the difference equations requires that the Courant 
condition be satisfied at each axial position X as, ^ 

“ - exp (-^'^) u o iR - °- 5 r 4R ) • (U3> ! 

x ' o 

I 

For the jet problem, the non-dimensional velocity u Q varies from one at 
the center to zero at the jet boundary. Therefore, it is difficult to satisfy 
the stability condition. This, however, can be ameliorated by specifying 
a small but non- zero value of u Q at the jet boundary. 

There exists another difficulty associated with the small u Q for 
R ■*> • . The difficulty arises because the gradient of the flux function, 

3F/3R , is small while the source function, S , does not vanish due to 
the terms containing l/u Q . This causes an Inaccuracy near the Jet boundary 
which then propagates Inward. This difficulty, of course, can be remedied 
by demanding a higher degree of accuracy for 0 by reducing the grid size, 

AR . However, the required computer time rapidly increases with the reduction 
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of AR . For Instance, halving Of the AR results in quadrupling of 

the computer time. Hie problem was resolved by the combination of reduction of 


AR to a practically feasible value and imposition of the condition, 

c£ * min £cj , (R -*■ 00 )J , (Al4) 

on the oxidant which diffuses from jet boundary to the core reg ; . With 
the above stability and accuracy considerations , solution of the governing 
equations is obtained by the two-step Lax-Vfendroff procedure. 

One additional point should be made. The dependent variables of Eqs. 

p 

(27a) - (30b) are c iQ ,. <Vc 1 > , and <V c^> . These are functions of 
z 1 , and Zg , as defined in Eqs. (33) - (35), which are, in turn, functions 
of z + and z“ via Eq. (37). It was found that the stability and accuracy 
criteria were best satisfied by employing c* and as the dependent 
variables of the equations rather than c io , <Vc i > , etc. 

Computer Program 

The computer program is written in PI/I language. A brief description 
of the procedures (subroutines in FORTRAN language) used in the code is given. 
The program structure is implicit so that numerous variables are implicitly 
passed through the procedures. 

1. Procedures 

JETFLAME — This is the main procedure which calls the necessary 
procedures and controls the step-size DX . The standard step-size DXJ3AVE 

4l 


is first tried, and if < 0 at any point in R the step-size is 
reduced by 80 percent. 

SOLVE — This procedure performs the solution by the two-step Lax- 
Wendroff numerical scheme. 

SOURCE — This is the procedure in which P(£2) and S(fl) are 
computed and the procedure SOLVE calls SOURCE. 

BCUNDRY — The boundary values for c* , c£ , c^ , c^ , c” , c“ , 
c d , and c" are specified at X = X Q . 

ERRFCTN — This procedure computes the error function with the 
polynomial fit, 

erf (x) = f a t n , (A15) 

n=l n 

where 


1 + 0.32759 X ’ 


(Al6) 


REACTION_RATE — The overall effective reaction rate. 


<u» R d R , 


is computed where 


<U> ‘ ^ f Ca ° Cb ° * + + erf ( V3 v o) 
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TRAFEZO — Integration is performed by a trapezoidal scheme. 

2. Inputs to the Code 

It consists of the procedure INPOT and a file INPUT 1. The variables 
for the input are as follows. 

E e * 1 - (uyup) 


ETAMAX 

maximum value of n 

XSTART 

X-value at which the initial values are specified 

XSTOP 

maximum value of X 

XOOT 

output control parameter; prints for every XOOT 

DR 

grid size in R 

E_HAT 

E 

MB_OVER_MA 

¥a 

MD_OVER_MA 


A1 

a l 

A 2 

a 2 

A3 

a 3 

A4 

a 4 

D 

y 

SIGMA 

a 


SPECIES ( I, J) 1-th species at J-th grid point where 



1 * 1 * ° a > 


2 * °a * 


3 * % 
4 -v c : 


5 

6 

7 

8 


c h 


'h ' 


(A19) 


Ihe variables listed— except SPECIES (# , 1) and SPECIES (#, MAXG) — are given 
in the procedure INPUTS. However, they are overridden and replaced by 

+ + -f -f 

the values in the file INPUT 1. The values of c , c^ , c d , and c h 
at the jet center and at R °° must be given in the file INPUT 1. 

3. Outputs from the Code 

Ihe final outputs are the average mass fractions, the turbulent 

A 

diffusion fluxes, and the local chemical reaction rate, <w> . 

A list of the other program symbols is given below. 


czo 

z o 

<V*cz> 

<Vz> 

<WP> 

•00 A 

J <u»R d R 
*o 

WP(R) 

<»>, (Eq. (A18)) 

ETA 

n 


4. 
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Program Listing 



e: 


<c 


on 

LU 

o 

O 

o 

<E 

QZ 

3 

Q. 

<3 



QZ 

<E 

LU 


H 


3 


CL 

S 

£3 

-J 

O 

CL 

o 

w 



* 
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FLAPIEJPROC OPTIONS(MAIN); 

DCL(SPECIES(*,*>, F(*,»), S(*,*), SPTC*,*)) CTl FLOAT DECU6); 
DCL (CHEM-REACTIONU), SPECIES-SAVEU,*) ) CTL FLOAT DECU6 >; 
DCL(AA, BB, Al, A2, A3, A4, 

X, X2, DX, DX-SAUE, XOUT, ETAMAX, 

XSTART, XSTOP, 

R, DR, RMAX, E, 

A, D, E-HAT, LAMDA, SIGMA, REF-CB0, 

MB-OVER-NA, MD-OUER-MA, SQRTE ) FLOAT DEC(IS); 

DCL (SET-PROFILE, VES INIT('Y'), NO INIT('N')) CHAR(1 ); 
DCL(GRID, DAXG, MAXR, INCREMENT, ITR, ITR-HAX ) FIXED BINC15); 
DCL (EXP, MAX, MIN, SORT, FLOAT, ABS, FIXED) BUILTIN; 

DCL( I, J ) FIXED BIN(15); 

DCL INPUT 1 FILE; 

/* SPECIES ( I )» C- UECTOR 
WHERE 1-CA0 

2- <U*CA>, 

3- CB0, 

4- <V*CB>, 

5- CD0, 

6»<U*CD>, 

?• H0, 

8-<U*H> */ 

CALL INPUTS; 

AA*SQRT ( 1 . S/E-HAT ) ; 

BB-SQRT ( E-HAT/18 . 848 ) ; 

DX-SAUE»0.?0*(1.-E)*DR; 

INCREMENT »F IXED ( 0 . 1/DR ) ; 

IF INCREMENT-0 THEN INCREMENT-1; 

X2,X»XSTART ; 

SQRTE-SQRT(E-HAT); 

CALL BOUNDRY( SPECIES); 

SPT-SPECIES; 

MAXR-MAXG-2; 


■ jsy. 


• V - • 


LAMDA*A4*<X+A2)/$IGnA; 1 

SET-PROF I LE ■ YES ; 

DX-DX-SAVE; 

MAXR»FIXED(20.*(X+A2)/(SIGf1A*DR ) ); 

CALL 0UT1; 

LLLJDO J-l TO 15; 

CALL SOLVE ( SPEC IES , F # S # RETRY ) ; 

END LLL; 

SET-PROFILE* NO; 

CALL OUT1; 

MAIN: DO 1<HILE<X<XSTOP ); 

DX-DX-SAUE; 

RETRY *X*X+DX; 

IF DX<1*E-09*DX-SAVE THEN 
DO; 

PUT SKIPC4) EDITC 'SPECIES< 1 ) IS NEGATIVE' ) 
<C0LC4),A); 

GO TO FINISH; 

END; 

LL2*D0 J* 10 TO NAXG; 

P1AXR* J; 

IF SPECIES<1/J)>0.9998*SPECIESC1/HAXG> THEN GO TO 0P2; 
END; 

0P2 : LAPJDA • A4* < X+A2 )/SIGPIA ; 

NAXR-P1AXR+8; 

CALL SOLVE ( SPECIES/ F/S # RETRY ) ; 

IF SPECIES <3/ 1 ) <0*005 THEN GO TO FINISH; 

IF X>*X2 THEN DO; 

X2-X2+X0UT; 

CALL 0UT1 ; 

END; 

IF SPECIESC3/1 ) <0.005 THEN GO TO FINISH; 

END WAIN; 


.. 








h 
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77. 

78* 

79, 

80. 
81. 
82. 

83. 

84. 

85. 
88 . 

87. 

88 . 

89. 

90. 

91. 

92. 

93. 

94. 

95. 
96 . 

97. 

98. 

99. 
100 . 
101 . 
102 . 

103. 

104. 

105. 

106. 

107. 

108. 

109. 

110 . 
111 . 






SOLVE : PROC ( SPECIES ,F,S, RETRY ) ; 

DCL(SPECIES(*,*),FU,*),S(*,*), SIGN, SLOPE) FLOAT DEC (16), 

(I, GRID, IX) FIXED BINC15); 

DCL J FIXED BIN<15); 

DCL RETRY LABEL; 

CALL SOURCE (SPECIES,F,S); 

FIRSTtDO GRID-2 TO NAXR-1; 

SPEC U DO 1-1 TO 8; 

IF ABS<FCI,GRID)»l.E-30 THEN 
SPT(I,GRID )-SPECIES( I,GRID)-(F(I,GRID+1 ) 

/F ( I, GRID )-l . X F(I,GRID XDX/DR 
♦DX*S(I,GRID> ; 

ELSE 

SPT( I, GRID ) -SPECIES (I, GRID) 

-(FU, GRID+1)-F(I,GRID))*DX/DR +DX*S(I,GRID) ; 

IF SET-PROFILE-NO THEN 
CHECKUDO; 

SPT(4, GRID )-NAX(0.0,SPT(4, GRID)); 

IF SPT(1, GRID) <0.0 THEN 
DO; 

X-X-DX; 

DX-0.2SDX; 

IF ABS(SPT(l,GRID))>l.E-06 THEN GO TO RETRY; 
ELSE SPT(l,GRID)-0.0; 

END; 

END CHECK 1 ; 

ELSE SPT(1, GRID )»HAXC0.0,SPTC1, GRID)); 

IF ABS(SPT(I,GRID) Xl.E-12 THEN SPTCI,GRID)-0.0; 

END SPEC1 ; 

END FIRST; 

DO 1-1 TO 4; 

SPT(2*I-1,1)«SPT(2*I-1,2); 

END; 

CALL SOURCE(SPT,F,S); 

SECOND t DO GRID-2 TO NAXR-1; 




112 . 

113. 

114. 

115. 

116. 

117. 

118. 
119. 
128. 
121 . 
122 . 

123. 

124. 

125. 

126. 

127. 

128. 
129. 
138. 

131. 

132. 

133. 

134. 

135. 

136. 

137. 

138. 

139. 
148. 

141. 

142. 

143. 

144. 

145. 

146. 


8 

8 

8 

8 



SPEC2SD0 I-l TO 8; 

IF ABS(F(I,GRID-1 ) )>1 .E-38 THEN 
SPECIES ( I , GRID > -8 . 5* ( SPT < I , GRID )+SPECIES< I, GRID ) 

-<F(I,GRID)/F< I,GRID-1 )-l. )*F< I, GRID-1 )*DX/DR 
+DX*S( INGRID)) % 

ELSE 

SPECIES(I,GRID)-8.S*(SPT(I, GRID >+SPECIES(I, GRID) 

-<F< I,GRID )-F< I, GRID-1 ) )*DX/DR 
+DX*S(I,GRID>) ; 

IF ABS(SPECIES(I,GRID ) Kl.E-12 THEN SPECIESU # GRIDl-8.8) 
END SPEC2; 

END SECOND) 

DO I»1 TO 4; 

SPECIES ( 211-1,1 )»SPECIES(2*I-i,2>) 

END) 

END SOL ME; 


SOURCE : PROC ( SPEC I ES , F , S ) ) 

DCL(SPECIES(*,* ), F<*,*), S<*,*>, 

ETA, U8, U0, U0.OUER.U8 , ERF , R, 

Gl, G2, G3, UA, 

F5, UC, 

CHEN1, CHEPI2, STCHHTRY, 

TERI11, TERN2, TERFI4, TERPI6, 

MSP (4 ), Ul, U2, 

DERI.M8, DERI.M8.X, DERI.U0.R,DERI-Ue.R,DERI.LOGUC, 
PHI_PLUS, PHI-NINUS ) FLOAT DECC16)) 

DCL SPECIU,*) CTL FLOAT DECC16)) 

DCL (I,J) FIXED BIN( 15 )) 

ALLOCATE SPECI(8,NAXR)) 

BOUNDRY t DO I-l TO 4) 

F < 2* 1 , 1 ) - 1 . 447282*SQRTE*SPECIES < 2*1 -1 , 1 )/2 . ) 

END BOUNDRY) 


147. 

148* 

149. 

150 * 

151. 

152. 

153. 

154. 

155. 

156. 

157. 

158. 

159. 

160. 
161. 
162. 

163. 

164. 

165. 

166. 

167. 

168. 

169. 

170. 

171. 

172. 

173. 

174. 

175. 

176. 

177. 

178. 

179. 

180. 
181. 


IF SET-PROFILE* YES THEN 
LLtDO 1*2 TO NAXR; 

IF SPECIESC1, 1X0.0 THEN SPECIE$(1,I )-0.0, 

IF SPECIES(2 # I )>0.0 THEN SPECIES(2,I >*SPECIESC2,I-1 ); 
END LLi 
L8SD0 1*1 TO 8; 

L6<D0 J-l TO HAXR; 

IF ABS(SPECIE$(I # J ) Xl.E-12 THEN SPECIESCI, J)-0.Oi 
END L6; 

END L8; 

R*0.0j 

Lit DO J«2 TO NAXR-lj 
R-R+DR; 

ETA*SIGNAXR/(X+A2); 

F5«l./(l.+0. 25XETAXETA ) j 

UC*Al/( A2+5.8 ) - EX(Ai/(A2+5.8) - A1/CA2+X) ), 
U0*1.-EX(1.-FSXX2); 

U0* EXA3X ( ETA-0 . 25XETAXX3 )XF5XX2; 

DERI-U0-R— EXETAXFSXX3XSIGMA/(X+A2>; 

DERI-V0*EXA3X< 1 .-0. 7SXETAXETA )XF5XF5-U0XETAXF5j 

DERI-V0.R-DERI.O0XSIGFIA/(X+A2)j 

DERI-U0-X— DERI-O0XETA/ CX+A2 ) ; 

DER I -LOGUC —EX A 1 / ( C X+ A2 ) XX2XUC ) ; 

VA-AAXU0; 

IF ABS(UA)>0.09 THEN CALL ERRFCTN(UA # 61); 

ELSE Gl-UA; 

G2-EXPC-VAXUA); 

G3-1./G2; 

CONTROL-AT-LARGE-ETA t 

IF ETA>2.45 THEN 

SPECIES<2,J>— ABSCSPECIES(2,J>); 

SPECIES ( 1, J )*MIN( SPECIES ( 1 # J ),SPECIES(1 # HAXG) ); 

L9tDO 1-1 TO 4; 

SPECI ( 2X1-1 # J ) *0 • 5X ( SPECIES ( 2X1-1 , J HI . 2247448/SQRTE 

XU0XSPECIES(2XI,J))j 


182* 

183. 

184. 

185. 

186. 
18?. 
188. 
189. 
198. 

191. 

192. 

193. 

194. 

195. 

196. 

197. 

198. 

199. 
280. 
281. 
282. 

283. 

284. 

285. 

286. 

287. 

288. 
289. 
218. 
211 . 
212 . 

213. 

214. 

215. 

216. 


SPEC I ( 2* I # J ) ■ 8 • 238 3294*SQRTE*G2*SPEC IES(2*I#J); 

END L9; 

L3*D0 I«1 TO 4; 

F(2*I-1,J)«2. * V8/U8 *$PECI < 2* I -1 , J )+2 . *SPEC I C 281 # J VU8 
>1 . ??3245385*V8*G3*SPEC I ( 2*1-1, J >/U8, 
F(2*I, J) - 1 . 4472825*SQRTE*G3*SPECI (2*1-1, J)/U8; 

END L3; 

IF SPECK1, J)-0.8 THEN SPECIC2, J)«8.8; 

IF SPECICl, JX5.8E-83 THEN 
DO; 

CHEN1 -SQRTE*D/LANDA*SPECI ( 1 , J ) , 

CHEH2-8.8; 

CHEfl_R£ACTION( J )-CHEf11 1 
GO TO OPT; 

END; 

ELSE IF SPECI ( 3, J ) <5 • 8E-83 THEN 
DO; 

CHEN1 -8 . 2S*SGRTE*D/LANDA*SPECI (3, J XNB-OVER-HA j 
CHEN2-8.8; 

CHEN-REACTIONC J J-CHEN1 ; 

GO TO OPT; 

END; 

ELSE L7JD0; 

TERN1 - ( 8 . 25* ( SPECIES ( 3, J )+ SPECIES ( A, J ) >**2 

+NB-OVER-NA* ( SPECIES ( 1 , J >+SPECIES (2,J)) 

*(SPECIES(3, J)4SPECIES(4,J) ) ) ; 
TERN2-(6.25*(SPECIES(3, J )-SPECIES(4, J) )**2 

♦NB-OVER-HA* C SPECIES C 1 , J ) -SPECIES < 2, J > > 

*(SPECIES(3, J)-SPECIES(4, J)) ) ; 

IF ABS < TERN 1 -TERH2 XI. £-86 THEN 
DO; 

PHI.PLUS-2. /TERN1; 

PHI.NINUS-8.8; 

END; 


""^.'Pippifw 


'■»tvwtrr' ^Tip^wp ' J 



21 ? • 
218 . 
219 . 
229 . 
221 . 
222 . 

223 . 

224 . 

225 . 

226 . 

227 . 

228 . 
229 . 
239 . 

231 . 

232 . 

233 . 

234 . 

235 . 

236 . 

237 . 

238 . 

239 . 
249 . 

241 . 

242 . 

243 . 

244 . 

245 . 

246 . 

247 . 

248 . 

249 . 
259 . 
251 . 


ELSE DO | 

PHI.PLUS- 1./TERH1 + 1./TERH2 % 

PHI-NINUS- 1 • /TERM - 1. /TERNS; 

END} 

CHEN1*HAX(9. , 9 . S9SQRTE9D/LANDA9SPECI < 1 , J >9SPECI (3, J >992 

9(PHI.PLU$+G19PHI.HINUS> > ; 

CHEN2-NAX( 9 .,9. 23932929E.HAT9D/LAHDA9SPECI (1,J J9SPECI C3 # J >992 

9G29PHI-HINUS ); 

CHEH.RE ACT I ON ( J ) -9 . 59 SQRTE9D/LANDA9SPEC I <i,J >9 SPEC I ( 3 # J)992 

9 ( PHI-PLUS+GltPHI-MNUS ) ; 

END L7| 

OPT * TERH4 - C U99DER I -U9-X 4- U99DERI-U9.R +V99U99DERI-L0GUC)t 
TERN6 - SGRTE/L AND AX ( 1 . *2.9D )/2.9j 
/* SOURCE TERNS FOR ORIGINAL EQUATIONS PSI'S 9/ 

S(1 # J>— 9.59CHEN1; 

S(2,J> — (TERN4*SPECI(1,J>+TERH69SPECI(2,JH8.59CHEII2>* 

S(3 # J )— NB-OUER-NAtCHENl j 

S(4,J> — (TERN4*SPECI(3,J>+TERN69SPECI(4,J) 

+NB-0UER.NA9CHEH2) ; 

S(5,J) -ND.OUER-N A9CHEN1 * 

S(6, J )--(TERN4XSPECI (5 # J )+SPECI (6# J )9TERH6 ) 

+ND.0UER.NA9CHEH2 ; 

S ( 7 ♦ J ) -HD.OUER JTIA9CHEH1 j 

S(8#J>— ( TERN49SPECI ( 7, J J4-SPECI < 8 # J J9TERN6 ) 

+ND.0UER-NA9CHEN2 ; 

L4SD0 1-1 TO 4} 

S(29I-1, J )• C-i .772453859(1 . 43 . 9UO9U0/E-HAT )*G39DERI_U9_X/U9 
+1 . 772453859U99G39DERI -U9-R/U0992 

4-2 • 9 ( DER I -U9-R/U0~DERI -U0-R9U0/U9992 ) ) 

XSPECK 2*1-1, J) ♦ (-5.317361S5/E-HAT9G3 
9 ( DER I -U9.X-U9XDER I -LOGUC > -2.9DERI-U9-R/U9992 
-2 . / < U99R ) -IS. 952984 6/E-HAT9929U99U99G39DERX-U9-X ) 
9SPECI (291, J ) 4>2.9S(29I-1,J)/U9 
-5 . 3173615539G39U99S ( 291 , J >/ ( E-HAT9U9 > | 
S(29I # J)-(4.34169?52/SQRTE9U9/U99DERI-U9-X 


252. 

253 . 

254 . 

255 . 

256 . 

257 . 

258 . 

259 . 
269. 
261 . 
262 . 

263 . 

264 . 

265 . 

266 . 

267 . 

268 . 
269 . 
279 . 

271 . 

272 . 

273 . 

274 . 

275 . 

276 . 

277 . 

278 . 

279 . 
289 . 
281 . 
282 . 

283 . 

284 . 

285 . 

286 . 


-1.4472925t$QRTEtG3XDERI.U9.R/U9X*2)8SPECX(2*I-l # J) 

♦ < -4 . 34169752XG3/SQRTEXDERI JLOGUC 
♦13 . 924822/ ( SQRTEXE.HAT )X09XG3SDERI_U9_X )1SPECI(28I,«I ) 

♦4 • 34169752/S0RTE8638S (2*1 , J »/U9 1 

END L4f 
END LI j 
FREE SPECI | 

END SOURCE} 

9 

9 

9 

9 BOUNDRV : PROC ( SPECIES ) j 

DCL( SPECIESU,*), ETA, R, TERN, U9.0UER.U9 ,U9> FLOAT DEC (16), 
DCLCEXPl, EXP2) FLOAT DEC(16); 

DCL (1, J) “IXED BIH(15), 

Ll*D0 , I-2 TO WAXG; 

R-R4DR} 

ETA*SIGHATR/(X+A2), 

TERPI • 1 . / ( 1 . +9 . 25XETAXETA )«2 , 

SPECIESC 1, 1 ) s SPECIES( 1,NAXG)$(1 .-TERf!)} 

SPECIES (3, 1 )-SPECIESC3, 1 > STERN} 

EXP1 - 1 . -EXP ( -9 . 999SSRSS6 ) , 

EXP2-EXPC-1.5XETA)} 

IF EXPK1.E-12 THEN EXP1-9.0} 

IF EXP2<1.E-12 THEN EXP2-0.9, 

IF ETA<13. THEN SPECIES(i,I )-SPECIES(i,I )XEXP1, 

IF ETA>7 .5 THEN SPECIES<3,I J-SPECIESC3, I )*EXP2, 

SPECIES (5, I )* (SPECIESC5, 1 )-SPECIES(S,HAXG) )STERN 

♦SPECIES (S, HAXG ) | 

SPECIES (7, I >■ (SPECIESC7, 1 )-SPECI£S(7 # HAXG) ) STERN 

♦SPECIES ( 7, NAXG )} 

END LI } 

F ( 4, 1 )- 1 . 4472925 SSQRTESSPECIES( 3, 1 )/2. } 
FC6,1>-1.447292SSSQRTEXSPECIESCS,1 >/2.j 


KjJ 


88 ?. 

888 . 

889. 

898. 

891. 

898. 

893. 

894. 

895. 

896. 
89?. 

898. 

899. 
380. 
301. 
308. 

303. 

304. 
36S. 

306. 

307. 

308. 

309. 

310. 

311. 
318 . 

313. 

314. 

315. 

316. 

317. 

318. 

319. 

380. 

381. 


0 

0 

0 

0 

0 


0 

0 


F<8,1 >-1.4472085*$QRTE*SPECIES(?,l )/8. • 
END BOUNDRY; 9 


ERRFCTN i PROC ( YY, ERF ) ; 

DCLCVV, V, ERF, P, T, EXPV8 ) FLOAT DECC16); 
DCL ACS) FLOAT DECC16) I NIT ( 0 • 254829598, -0 


DCL I FIXED BIN(IS); 
Y-ABS(VY); 

IF V<0.01 THEN 


1.481413741, -1 
1.061405489 ) ; 


.884496736, 

•453152027, 


DO; 

ERF-0.0; 

RETURN; 

END; 

IF Y>4.0 THEN DO; 

ERF-1.0; 

GO TO RTN; 

END; 

P-0.32759; 

T-l./Cl.+PJY); 

ERF-0.0; 

DO I-l TO 5; 

ERF-ERF+ACI )$Tt*I; 

END; 

EXPY8-EXPC-YtY>; 

IF EXPY8<1.0E-08 THEN EXPY8-0.0; 
ERF-1.-ERFTEXPY2; 

RTN: IF VY<0.0 THEN ERF— ERF; 

END ERRFCTN; 


322 . e 

323. 0 REACT ION-RATE* PROC (UP); 

324. DCLOJC, UP, RESULT, ETA,R,F5,FCTN(* ) CTL) FLOAT DEC<16); 

325. DCL J FIXED BIN(15); 

326. ALLOCATE FCTNOIAXR); 

327. R-O.0> 

328. FCTN( 1 )-0.O; 

329. LI* DO J-2 TO HAXR-1; 

330. R-R+DR; 

331. UC-Al/( A2+5.8 ) - E*( Al/( A2+5.8 ) - A1/(A2+X) >; 

332 . FCTN < J ) -UCtCHEH-REACTION ( J >*R; 

333. END LI; 

334. CALL TRAPEZO(FCTN,DR,HAXR, RESULT); 

335. UP -RESULT; 

336. FREE FCTN; 

337. END REACTION-RATE; 

338. 0 

339 0 

340 ! 0 TR APEZO * PROC ( FCTN , DR , MAXR , RESULT ) ; 

341. DCL(FCTNU), DR, RESULT) FLOAT DECC16); 

342. DCL< I,I1AXR ) FIXED BIN(iS); 

343. RESULT-0.0; 

344. LL*DO 1-1 TO PIAXR-2; 

345. RESULT«RESULT+(FCTN< I )+FCTN( 1+1 ) >*DR; 

346. END LL; 

347. END TRAPEZO; 

348. 0 

349. 0 

350. 0 

351. 0 INPUTS: PROC; 

352. DCL AAAA CHAR(l), BBBB CHAR(l); 

353. E-0.80; 

354. ETANAX-8.0; 

C55. XSTART-5.8; 

c56. XSTOP-100. ; 



35?. DR-0.20; 

358. D-0. 333333333; 

359. E-HAT *0.10; 

369. MB.OUER .MA-2 • /32 • ; 

361. ND-OUER-HA* 1 8 . /32 . ; 

362. Al-6.4; 

363. A2-6.6; 

364. A3-0. 03295; 

365. A4-0.10; 

366. SIGMA-16.; 

367. GET DATAt 

368. E 



? 



I' 




f 

!■ 


* 

li 

f; 


369. 

370. 

371. 

372. 

373. 

374. 

375. 

376. 

377. 

378. 

379. 

380. 

381. 

382. 

383. 

384. 

385. 

386. 

387. 

388. 

389. 

390. 

391. 


.ETAMAX 

.XSTART 

.XSTOP 

,DR 

.XOUT 

.E-HAT 

.MB-OVER-HA 

.MD-OUER-MA 

A4 

) FILE(INPUTl) COPV(SVSPRINT); 

MAXG-F IXED ( ETAMAX* ( XST0P+A2 ) / ( SIGMA*DR ) ) ; 

ALLOCATE SPECIES(8,MAXG). 

SPT(8.P1AXG )> 

F<8,MAXG), 

CHEM_REACTION(MAXG ). 

S(8.MAXG); 

SPECIES. SPT . F. S-0.0; CHEM.REACTION -0.0; 

GET LISTCAAAA. BBBB 

.SPECIES(l.l), SPECIES (l.HAXG) 

. SPECIES(3. 1 ). SPECIES ( 3, MAXG) 

.SPECIESC5.1 ). SPECIES ( 5. MAXG) 

.SPECIESC7.1 ). SPECIESC7.HAX6) 

) FILE( INPUT1 ) COPY ( SVSPRINT ) ; 




392. 

393. 

394. 

395. 

396. 

397. 

398. 

399. 

400. 

401. 

402. 

403. 

404. 

405. 

406. 

407. 

408. 

409. 

410. 

411. 

412. 

413. 

414. 

415. 

416. 

417. 

418. 

419. 

420. 

421. 

422. 

423. 

424. 

425. 

426. 




END INPUTS; 

0 
0 
0 

0 0UT1JPR0C; 

DCL ( R, UP. UPR, UC, CD, ETA, F5, U0, 62) FLOAT DECC16); 

DCL SPECK*,*) CTL FLOAT DEC (16); 

DCL <J,I) FIXED BIN< 15 ); 

ALLOCATE SPECI (8,J1AXR ); 

PUT PAGE FILE( SYSPRINT ) EDIT('X»',X) (C0L(3),A,C0L(6),F(9,5>); 
CALL REACTION_RATE(UP ); 

PUT SKIP FILE(SYSPRINT) EDIT( * <UP>* / ,UP)(COL(3),A,E(10,3)); 

PUT FILE(SYC?RINT ) EDIK 'R', 'CA0', '<U*CA>', 'CB0', , <U*CB> / , 
'CD0', '<U*CD>', ' H0', '<V*H>','CE0', 'UP(R)', 'ETA' ) 

(COLO), A, COLO ), A, C0L(19),A, C0L(29),A, C0L(39),A, 
C0L(49 ), A, C0L(59 ), A, C0L(69),A, C0LC79 ), A,C0L(91 ),A, 
COL( 102), A, COL ( 113), A); 

R*0 # 0j 

DO 1-1 TO MAXR BY INCREMENT ; 

ETA® SIGN A*R/'(X+A2 ) ; 

F5«1./(1.+0.25*ETA*ETA); 

V0» E*A3* ( ETA-0 . 25*ETA**3 )*F5**2 ; 

G2®EXP<-1.5*U0**2/E-HAT); 

LL3JD0 J-l TO 4; 

SPEC I ( 2* J- 1 , I ) *0 . 5* ( SPEC I ES ( 2* J-l , I ) 

+ 1 . 2247448/SQRTE*U0*SPECIES(2*J, I ) ); 
SPECI ( 2* J, I ) -0 . 2303294*SQRTE*G2*SPECIES ( 2* J, I ) ; 

END LL3; 

CD-PIAX(0.0, 1 . -SPECI <1,1 )-SPECI (3, 1 )-SPECI (5, 1 ) ); 

UC-Al/(A2+5.8) - E*< Ai/< A2+5.8 ) - A1/(A2+X) ); 
UPR-UC*CHEN_REACTION< I ) ; 

PUT FILE(SYSPRINT) EDITCR, SPECI(1,I), SPECI(2,I), 

SPECI (3,1), SPECI (4, 1 ), SPECKS, I), SPECI(6,I>, 

SPECK7, 1 ), SPECK 8, 1 ), CD, UPR, ETA ) 

<COL( 1 ),F<5,2 ), C0L(7),E(9,2), C0L(17),EC9,2), C0L<27),EC9,2>, 




427. 

428. 

429. 

430. 

431. 

432. 

433. 

434. 

435. 

436. 


hiwiww,! .1 


IN qpp 




C0L<37),E<9,2), C0l<47),E<9,2>, C0LCS7 ),E<9,2 ), C0L(67),E<9,2), 
COL(77),E(9,2)*COL(88) # E(9,2),COL(99),E(9 # 2),COL(110),F(7 # 3) ); 
R«R*fINCREHENT*DR; 

END; 

FREE SPEC I; 

END 0UT1; 

0 

0 

0 FINISHtCALL 0UT1; 

END FLANE; 



//G0.INPUT1 DD * 
E-0.80 
ETANAX-8.0 
XSTOP-60. 
DR-0.20 
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Axisynmetric turbulent jet 









Sketch of characteristics 
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Figure 3.- Initial profiles specified. 
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Figure 13. Axial variations of combustion rate. 




Figure 14.- Axial variations of corrtoustion rate. 
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Figure 15.- Center line variations of mean mass fractions and tenperature. 
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Figure 16.- Center line variations of mean mass fractions and 
temperature . 
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^Flame Sheet Sol'n 
E“0.1, a4”0.5 

u<x)/up“0.2 
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Figure. 17.- Center line variations of mean mass fractions and terrperature 











Profiles at X =5.8 
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^Flame Sheet Sol’n 
E = 0 . 1 , 84=0.5 
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